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Abstract 

Evolutionary branching is analysed in a stochastic, individual-based population model under mutation and selection. In such 
C^l models, the common assumption is that individual reproduction and life career are characterised by values of a trait, and also by 
population sizes, and that mutations lead to small changes e in trait value. Then, traditionally, the evolutionary dynamics is studied 



in the limit e — > 0. In the present approach, small but non-negligible mutational steps are considered. By means of theoretical 



analysis in the limit of infinitely large populations, as well as computer simulations, we demonstrate how discrete mutational steps 
CJ affect the patterns of evolutionary branching. We also argue that the average time to the first branching depends in a sensitive way 
on both mutational step size and population size. 

1— I Keywords: Evolutionary branching, genetic drift, selection, adaptation 
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1. Introduction 

The development of populations subject to frequency- or 
density-dependent selection is an important and central but not 
easily analysed topic in theoretical biology. Many different 
models and approximation schemes have been developed in or- 
der to understand how trait values change and how the intrigu- 
ing possibility of 'evolutionary branching' may arise, where 
a population of individuals with a well-defined trait value or 
genotype (monomorphic, in other words) experiences a split 
leading to two or more coexisting trait values or genotypes in 
the population. This process is important because it gives rise 
to increased biodiversity. 

A variety of population-genetic and game-theoretic meth- 
ods, as well as the body of ideas and approx i matio n 
schemes kn o wn as a daptive dynamic s, cf. Metz et al. J 1996b; 
Gerit z et al. dl998h : Diekman (2004); Waxman and Gavrilets 



X 



'"J d2005h : IMetz et all d2012h . have been employed in order to 



study evolutionary branching. Several simplifications are com- 
monly resorted to. First, complications due to mating and re- 
combination in diploid populations are disregarded by ascrib- 
ing to each individual a single trait value x that is transferred by 
clonal reproduction unless mutation occurs (see however the re- 
cent interesting approach to adapti ve dynamics with Mendelian 



inheritance by ICollet et al.L 1201 ll) . Selection is taken to act 



through a trait- and density-dependent fitness function, identi- 
fied with mean reproduction. Second, it is assumed that fitness 
is a smooth function of the trait value x and the population den- 
sity. The third postulate is that the mutation rate // is so small 
that only few trait values are represented in the population at 
any one time. Mutations thus occur so rarely that we can talk 



of a separation of time scales between long-term evolutionary 
dynamics (the sequence of mutations) and short-term popula- 
tion dynamics, defined by birth and death events, an evolution- 
ary versus an ecological time scale. Fourth, it is assumed that 
mutations lead only to small changes in the trait value (that mu- 
tational step sizes e are small). Fifth, the population size N is 
taken to be large, and effects of random genetic drift are ne- 
glected (apart from incorporating the fixation probability of an 
advantageous mutation). 

In the adaptive-dynamics literature a further step is taken 
by letting N — » oo, /u — > 0, and e — > 0. Then evolution 
of a monomorphic population to the first branching of traits 
is analytically described by the so-called 'canonical equation' 



( Dieckmann and Law . 1996), or more generally by determin- 



istic evolution on the fitness landscape dWaxman and GavriletsL 



2005al : lLande , llwasa et alll 199 11) . How these determinis 
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tic adaptive dynamics may arise from a strict stochastic popula- 
ti on model has been de scribed in detail in the elegant treatment 
of lChampagnal d2006h and later papers 

As a monomorphic (single-trait) population approaches a lo- 
cal fitness maximum where evolutionary branching may occur 
in the sense that populations with two different trait values may 
coexist from there onwards, the canonical equation fails to de- 
scribe the dynamics. It becomes necessary to ask how stochas- 
ticity in a finite population affects the possibility of evolutionary 
branching. Further questions are: what is the shape of the bifur- 
cation diagram of trait values in the wake of a branching? How 
do the values of ft and e influence the evolutionary dynamics 
in a finite population? How long does it take (in the mean) for 
evolutionary branching to occur? 

We address these matters analytically and by computer simu- 
lations of a stochastic, individual-based model for the evolution 
of trait values in a finite population subject to density-dependent 
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Figure 1 : A schematic view of the deterministic model. The left panel depicts 
the first three evolutionary branching events. The right panel illustrates the 
main ingredients of the fitness function given by Eq. j5J. The function is 
plotted for a particular value of y = xq which in the case shown is positive. 



fitness function symmetric with respect to its maximal value at 
the trait x — 0. In the limit of N — > oo, p — > 0, and e small, we 
show that the evolution of a monomorphic population towards 
the fitness optimum at x — is governed by a canonical equa- 
tion. For our model, standard stability analysis in the limit of 
e — > would predict evolutionary branching of a monomorphic 
population at x — 0. We demonstrate that at small but fixed 
mutational step sizes e, evolutionary branching typically occurs 
at non-zero values of x. In the limit of infinite population size 
and vanishing mutation rate we compute the critical value of x 
where the first branching is most likely to occur, as a function 
of e. In addition we show that the evolutionary dynamics in 
this limit gives rise to a non-symmetric bifurcation diagram af- 
ter the first branching (for a symmetric model such as ours, the 
standard theory predicts a symmetric bifurcation diagram in the 
limit e — > 0). Moreover, we state geometrical conditions deter- 
mining the trait values where the second evolutionary branching 
event is most likely to occur. 

Finding the average time until evolutionary branching occurs 
in our model requires a more refined analysis. It is necessary to 
specify how N ^ oo and e — > 0. This determines the stochas- 
tic population dynamics in the critical stage preceding the first 
branching. That stage must last long enough to render branch- 
ing possible. We discuss results of computer simulations relat- 
ing to this question briefly, but a precise mathematical descrip- 
tion of the problem is left for future work. 

The remainder of this paper is organised as follows. In Sec- 
tion[2]we introduce the stochastic, individual-based model used 
in simulations. Section [3] summarises standard results obtained 
in the limit N — > oo, ^ — > 0, and e — > 0. These are con- 
trasted with numerical experiments at relevant values of N, /u, 
and e in Section [4] Discrepancies between the e — > predic- 
tions and the results of numerical experiments are discussed. 
Our results (and the calculations supporting them) addressing 
the issues raised in Section|4]are summarised in Section|5] Sec- 
tion [6] contains our conclusions. Three appendices summarise 
details of our calculations. 



2. Model 

The stochastic, individual-based model used in the numeri- 
cal experiments discussed below describes a finite population 
of individuals competing for resources. Each individual is as- 
signed a trait value x. Individuals with the same trait value (be- 
longing to the same ecological 'niche') compete for a limited 
amount of resources. The scarcity of resources is parametrised 
by a trait-dependent carrying capacity K x . There is an 'optimal' 
trait value corresponding to the value of x where K x assumes 
its global maximum. Individuals belonging to different niches 
(corresponding to different trait values, x and y, say) may also 
interact. The strength of interaction depends upon the differ- 
ence |jc - y\ in trait values. 

Mutations give rise to changes in the trait value, x — * x ± e. 
Here e > is a fixed mutational step size and the set of possible 
traits occurring in the population is x — 0, ±e, ±2e, Individ- 
ual fitness is measured by the mean offspring number, which is 
determined by the trait values and the intensity of competition. 
The latter, in its turn, is controlled by the scarcity of resources 
and the number of individuals. In other words, the model de- 
scribes the evolution of a trait value, or of trait values, in a pop- 
ulation subject to density-dependent selection. The population 
would be expected to evolve towards the optimal trait value. In- 
tense competition between individuals with trait values close to 
the optimal one may however give rise to sub-populations with 
well-defined but distinct trait values. This phenomenon is re- 
ferred to as evolutionary branching, and is the subject of this 
paper. 

The number of individuals with a given trait value x is de- 
noted by Z x . It is as sumed that each indiv idual has either 
none or two offspring (IKlebaner et all 1201 II) . The offspring 



of individuals in one generation constitute the next generation. 
We denote the fitness function by M t (Z v , y e D) which is the 
mean offspring number for an individual with trait value x, and 
where D is the set of trait values represented in the population. 
This fitness is a function of the co-existing sub-population sizes 
(Zy,y 6 D). We take it to be of the form: 



M x (Z y ,yeD) = 



l+(NK x )-i£ v£D y xy Z y 



(1) 



Here the carrying capacity of the subpopulation with trait value 
x is NK X , where N is a population-size parameter to be thought 
of as large. Competition between subpopulations correspond- 
ing to different trait values x and y is regulated by the function 

Jxy 

The offspring usually inherit the parental trait x, but with a 
small probability jj they mutate (independently) either to x— e or 
to x+ e, with equal probabilities. Since time is measured in gen- 
erations, the model is thus discrete both in time and trait space. 
It can b e viewed as the simpl est possible structure ('the bare 



bones', IKlebaner et al.L 1201 11) behind adaptive dynamics, and 



as such is an interesting mathematical object in its own right. 
Since all individuals have the same life span, one, it displays 
an unyielding age dependence which can be seen as a diametri- 
cal counterpart and natural complement to the equally artificial 
non-aging individuals of birth-and-death processes, underlying 



2 



the sequence of papers bv lChampagnat and Meleardl(l201 II) and 
others, and also classical determinist ic approaches. For more 
realistical ly age-structured models cf. Meleard and Tran ( 20091) 
as well as Lfagers and Klebanei ( 201 1 ). 

The form (Q} is certainly ad hoc, and should be viewed as an 
explor ative illustration. But it is not wi thout historical prece- 
dence (IChristianse n and Loes chke. Il980h . As there, we take K x 
and y n , to be Gaussian functions. Indeed, this is a choice often 
encountered in the adaptive-dynamics literature, 



trait to establish itself in the set of coexisting traits. If M? < 1, 
the mutation can not establish itself and is wiped out by genetic 
drift. 

In summary, we have chosen an extremely simple model 
which however retains so much of the structure of evolution- 
ary dynamics that the questions raised in the introduction can 
still be formulated. It is described by only four parameters, 
the population-size parameter N, the mutation rate fi, the muta- 
tional step size e, and the interaction parameter a. 



K x = e 



and 



7xy = e 



-(x-y) 2 (l+a) 



(2) 



The parameter a regulates the competition between subpopula- 
tions with different traits. A larger value of a implies weaker 
competition. We take a > 0. This ensures that competition be- 
tween subpopulations is weak enough to allow for evolutionary 
branchings. It is convenient to replace population sizes Z x by 
densities f x = Z X /(NK X ), so that the fitness function ([1]) takes 
the form 



M x (f y ,yeD) 



1 + TjyeD Rxyfy 



Here 



R xy = y xy K y IK x = e 



-(x-y)(ax-(2+a)y) 



(3) 



(4) 



is the fitness kernel which combines the effects of the interac- 
tion jxy and the carrying capacity K x . 

Eqs. © and © allow for explicit formulae for the equilib- 
rium densities (fi D ,y e D) at which subpopulations may co- 
exist in the absence of mutations. Note however that even 
this coexistence is temporary, al beit potentially of long dura- 
tion Pagers and Klebanen, 1201 lb . We shall still refer to them 
as equilibrium densities. They are defined by the condition 
M x = 1, where 



M D = 



D ' 



(5) 



The equation = 1 describes the critical regime of reproduc- 
tion. It corresponds to the condition 



yeD 



1, x e D. 



(6) 



In the monomorphic case D = {x} we find that M x (f x ) — 2/(1 + 
f x ) and in the absence of mutations, the density f x stabilises at 
the value f x = 1 . Substituting this result into Eq. (0 we obtain 



Ml xi 



I + e -(z-x)(az-(2+a)x) 



(7) 



Importantly, the function M? given by Eq. (f5J) plays a crucial 
role in the invasion analysis. Suppose an individual of type 
z i D appears in a stable population with the set of traits D as a 
mutant of x e D. Its fitness has the form 



M z {f?,yeD-j:) = 



Mt 



because R u = 1 and the small initial frequency of the mutant 
f* can be neglected. If M D > 1, there is a chance for the new 



3. Standard predictions in the limit e -» 

In spite of its simplicity, the model introduced in the previ- 
ous Section is not easily analysed. As pointed out, t he stan- 
dard adaptive-dynamics approach ( Geritz et all 1998 ) studies 
continuous-time population models in the limit Af — > oo and 
fx — * 0, and then e — > 0. In this Section we state the corre- 
sponding results for our discrete-time model. 

If can be thought of as infinite, law-of-large-number ef- 
fects replace random processes by their expectations, and if 
both jj — > and — > oo, then only a small number of trait 
values are represented in the population at any specific time: 
the basic cases studied below are those of monomorphic popu- 
lations, D = {x}, and dimorphic ones, D = {x\,X2}. 

In this limit, a time-scale separation appears between the evo- 
lutionary dynamics (referring to long-term changes of the trait 
value) and the ecological time scale of population dynamics. 
Letting e — > makes it possible to describe the dynamics of 
the trait value x by local stability analysis, central to adaptive 
dynamics. 

Fig. [TJ illustrates the standard predictions obtained in these 
limits. First, starting from an initially positive trait value, x(f) 
is expected to evolve deterministically towards the local op- 
timum (at x — in our model). The path x(t) is described 
by a differential equation referred to as a 'canonical e quation' 
(IDieckmann and Law . 119961: IChampagnat et all 1200 lb and re- 
lated to similar equation s studied in population genetics (Lande, 
llwasa et all 1 199 II) . Second, the first evolutionary branch- 
ing occurs at x = 0. Third, after the first branching, the pair of 
traits evolves deterministically in a symmetric fashion, and the 
second and third branchings occur simultaneously. 

3.1. Deterministic evolution towards the first branching 

Assume that the population starts from a positive initial trait 
value xq = joe which is large compared to the small value of 
the mutational step size e: jo » 1 . Classical adaptive dynamics 
then predicts that the initial trait value is consecutively replaced 
by smaller values of x, driving x(t) towards x = 0. The resulting 
function of time x{t) is described by a differential equation for 
x(t), known as the canon ical equation (of trait evolut ion), the 
standard r eference being Die ckmann and Lawl (119961) . In the 
notation of lChampagnat et al. I d200ll see Eq. (7a) ) the classical 
canonical equation of trait evolution has the form 



-jr = V(x)^-n(x)dif(x, x) . 
at 2 



(8) 



In this expression, fi(x) is the mutation rate for the trait value 
x, <x 2 is the mutational variance, n(x) is the equilibrium popula- 
tion size, and f(z, x) is the fitness function. Recall that within a 
model formulated in terms of stochastic birth and death pro- 
cesses, the fitness function f(z,x) is the difference between 
the birth and death rates for rare mutants with trait z in an x- 
population at its quasi-stable size. 

In our case, the mutation rate fi(x) = ju is independent of 
the trait value x, <r\ is simply e 2 , and n(x) = Ne . The fit- 
ness function f(z, x) corresponding to our discrete-time model 
should be calculated as 



/(z,*) = logM w , 



(9) 



where Mf* is given by Eq. (0. This formula relies on a well- 
known property of a linear birth-death process stemming from 
a single individual: if / is the difference between the birth and 
death rates per individual, the population size at time t has mean 
e'l . Eq. (O stipulates that the expected population size e^ z ' x ^ 
after one unit of continuous time is equal to the expected size 
Mf' of the first generation in our discrete-time model. 

From Eqs. (|9]i and © we conclude that the fitness gradient is 
given by 

dif(x,x) = -x, 

since Mi = 1. In our case we would thus expect the canonical 
equation to take the form 



dx e 2 . 
— = -u—Ne 
dt 2 



(10) 



However, the correct canonical equation for our model Eq. (TT~9b 
derived in Section IBTTl reveals that the trait substitution process 
goes twice as fast in the discrete-time model compared to the 
continuous-time model. This is reminiscent of the well-known 
phenomenon of genetic drift running twice as fast in the Moran 
model th an in the Wright -Fisher model of the same size, see for 
example IWakelevI (120081) . though the latter phenomenon has a 
different explanation. 

3.2. Location of the first evolutionary branching 

The canonical equation suggests that the substitution process 
of the trait value slows down as the point x — is approached. 
In the limit of N — » oo, ^ — » 0, and e — > 0, standard adaptive- 
dynamics analysis predicts the first evolutionary branching to 
occur at x — 0, giving rise to two branches, x\ < and X2 > 0. 

3.3. Second and third evolutionary branchings 

Moreover, in the very same limit, the deterministic evolution 
after the first branching is predicted to be symmetric. In other 
words, the two trait values x\ and X2 after the first branching 
are expected to evolve as x\{t) - —Xzif). The second and third 
evolutionary branchings are expected to occur simultaneously 
when X2{t) reaches a critical value, which in our case can be 
computed explicitly as 



log(2a+ 1) 



(ID 



1 + a 

As a function of a the value of x a reaches its maximum 
0.3731... ator= 1.2955.... 




Figure 2: a Shows ten realisations of the evolution of x(t) from an initially large 
value xq = 2 for a mutational step size of e = 1CT 2 (additional parameters: 
N = 10 6 , n = lfT 8 , a = 9, = fiN = 1CT 2 ). b Shows the trait value lasa 
function of the average time of staying at x, conditional on that no branching 
has occurred (dots). Also shown is an approximate solution of )10t (solid line). 



4. Computer simulations 

Figs. |21 [3]and|4]summarise results of direct numerical simu- 
lations of the model described in Section [2] for large values of 
N, small values of fi, and for a small value of e (equal to 10~ 2 ). 
In the algorithm, binomial distributions of numbers of mutants 
have been approximated by the appropriate Poisson distribu- 
tions. 

We now compare simulation results with the predictions from 
Section[3] Fig.|2^ shows ten realisations of the evolution of x(t) 
for the mutational step size e = 10~ 2 from an initially large 
value xq = 2. (The additional parameters are given in the figure 
caption). We see that x{t) decreases towards the first branching 
(which in all cases occurs for positive values of x and not at 
x = 0). Fig. |2]b shows a plot of the trait value x versus the 
average time of staying at x, conditional on no branching having 
occurred. Also shown is an approximate solution of the correct 
canonical equation derived in Section 15.11 as Eq. (fT9l l. This is 
expressed in terms of the expected time for x(t) to reach level 



t(x) = J] A; 1 



(12) 



where A x = fiNexe~ x ~ (see below) is the rate of the asymptot- 
ically exponential holding time at level x. For large values of 
x we observe good agreement between (fT2l and the average of 
ten independent realisations of the substitution process. But as 
x — is approached, the numerical average falls below the pre- 
dicted average. The neighbourhood of the first branching event 
is not described by the canonical equation. 
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Figure 3 : Results of computer simulations of trait evolution in the the stochastic 
individual-based model described in Section [5] Parameters: a = 9, e = 0.01, 
and a N = 5 x 10 6 , n = 2 x lO" 10 . 

Fig. [3] shows results of two computer simulations starting in 
the vicinity of x = 0. In Fig. [3^ the first branching does not 
occur at x = but at x = 0.05. Further, evolution after the 
first branching is not symmetrical. This can also be seen from 
Fig. 3] which displays the data of Fig. [3^ in the X1-X2 plane (the 
standard adaptive-dynamics approach predicts, as pointed out 
above, that the pair (xi(f), X2(t)) moves on the diagonal shown 
in Fig. |4] as a dashed line). Note also that the second and third 
branchings do not occur simultaneously, and not symmetrically 
in the simulations. 

Fig. 0> shows results of a direct numerical simulation with 
the same overall substitution rate 9 = fiN as Fig. [3^, but for 
a smaller value of N and a larger value of //. In this example, 
branching did not occur during the simulation time 6t < 10 5 . 
By contrast, the trait value is seen to diffuse around x — 0. 

In summary we find that the standard adaptive-dynamics ap- 
proach describes the initial time evolution of the trait value x{t) 
well, provided x(t) is large enough. In the vicinity of and after 
the first branching, however, the observed patterns of evolution- 
ary dynamics (Figs. [3]and|4]i differ from the standard predic- 
tions. In the following we show that the asymmetric patterns 
observed in Figs. [3j» and [4] are due to the positive mutational 
step size e. In the limit N — » 00 and fi — > we were able to 
characterise the branching patterns at small but fixed values of 
e. To determine how long it takes on average until evolutionary 
branching occurs in a finite population requires a more refined 
analysis. We return to this in Section IB~3~1 

5. Results 

5.1. Deterministic evolution towards the first branching 
In the monomorphic case (D = {x}) Eq. (0 implies 




-0.6 -0.4 -0.2 

Xi(t) 

Figure 4: Shows the trajectory of Fig. [3^ after the first branching in the X1-X2- 
plane. Only points corresponding to X\ < are shown. Parameters the same as 
Fig.[2». 



for z ~ x. For x > it follows that 

> 1, M^_ e < 1, M l *~ €] < 1, and M x x+e] > 1 . (14) 

These inequalities show that the expected trend of adaptive evo- 
lution in the monomorphic phase is simple: the initial positive 
trait value is consecutively replaced by smaller values driving 
x(t) towards zero. The dynamics of this trait substitution pro- 
cess is described by the canonical equation whose correct form 
for our model is derived next. 

Observe that the expected number of mutations x — > x - 
e in the x-population, assumed to be at its quasi-stable size, 
Z x - NK X — Ne , is (fi/2)Ne~ x per generation, where \x is 
divided by 2 since we discard the x — > x+e-mutations. Thus the 
expected waiting time between two consecutive replacements is 



^Ne- 
2 



(15) 



Here P : (x) stands for the probability that a population from one 
z-individual survives and takes over from the monomorphic x- 
population. To find the probability P-(x) of the mutant z = 
x - £ invading and replacing a resident population with trait 
x we can approximate the population dynamics of the invader 
by a binary branching process starting from a single individual. 
The survival probability of this branching process is 



P z (x) = 2- 



Ml 



1 



Ml' 



(16) 



1 * x(x-z) + (a/2-x 2 )(z-x) 2 



(13) 



see for example Eq. (5.66) in lHaccou et al.l ( 12005b . Eq. ( [TBI 
implies 

M£* e *l+ec, (17) 

and we obtain P x - e (x) ~ 2ex, at least for x far from zero. (The 
corresponding classical result for the Wright-Fisher model is 
the following. Suppose that an advantageous allele is intro- 
duced in the population with its mean offspring number being 
larger compared to the wild type by a factor ( 1 + s). Then the fix- 
ation probability of this allele is approximately given by P ~ 2s 
for small values of s.) 



Substituting P x - e (x) w 2ex into Eq. (TT3T > results in A x w 
fieNe'^x. Now, since the evolution of the trait value x(t) as 
a function of time is approximately given by 



df t(x - e) - t(x) 1 
dx e eA x ' 

we conclude that the canonical equation takes the form 



Ax 



— = -ue^Ne 
dt 



(18) 



(19) 



which is similar (but not identical) to Eq. dTOb . Indeed, Eqs. ( [19b 
and ( [Tol l differ by a factor of I. 

This factor c an be explained as follows. The standard deriva- 
tion o f Eq. © ( Dieckmann and Law , 19961 : Champagnat et al 



20011) refers to birth-and-death processes in continuous time (or 
corresponding deterministic formulations), whereas our model 
considers binary splitting in discrete, non-overlapping gener- 
ations. Our form of the equation appears expli citly for Pois- 
son reproduction in discrete time in iMetzl (1201 lb . The fact that 
different effective populati on sizes can appear in the canonica l 
equation was first noted bv lDurinx. Metz. and Meszenal (120081) . 

Note that the speed of the trait-substitution process given by 
Eq. © would not change if we were to take different values of 
the birth and death rates, as long as their difference f(z, x) re- 
mains the same. We argue that the proper choice of the birth 
rate in the monomorphic stable population is unity, matching 
the birth rate in our discrete-time model: one birth per gen- 
eration per individual, Mjf' = 1. Then, in accordance with 
Eq. (01, the corresponding survival probability of a single mu- 
tant is given by f(x - e, x) = log M^l e ~ ex. Comparing this 
with its discrete-time counterpart P A _ e (x) » 2ex, we conclude 
that the survival probability in the discrete model is two times 
larger. 

An important assumption behind Eq. ( fT9] i is that successful 
mutations should not come too soon after each other in order 
to make sure th at selective sweeps do not o verlap ('clonal in- 



terference'), see Su-Chan Park et al. ( 2010t) . As a very rough 



approximation, we estimate the average fixation time from the 
equation 



N. 



(20) 



This equation simply suggests to neglect the competition 
among individuals and equate the target population size N with 
the mutant reproduction factor for Tr x consecutive generations. 
Combining Eqs. $1% and (1201 we find T^ x « log N/ (ex) in the 
absence of competing mutations. The average time between 
successful mutations is A" 1 . So the condition ensuring non- 
overlapping sweeps is that the product A x T^x w fiN log Ne~ x ~ 
must be very small so that 



fiNlogN «: 1 



(21) 



In other words, condition Eq. (EH guarantees that after a mu- 
tation the new set of types settles in an equilibrium before the 
next successful mutation event. 
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Figure 5: Shows histograms of the locations j of the first branching event. All 
histograms are normalised to unity. Initial condition j = 7. Parameters: a = 9, 
e = l()~ 2 . The first evolutionary branching is expected to occur in the interval 
— 5 < j < 5. The upper limit of this interval is indicated by black arrows. 
Panel a: N = 5 X 10 6 , fi = 2 X lfT 10 , all 100 runs resulted in an evolutionary 
branching. Panel b: N = 10 6 , fi = 10" 9 , out of total 100 total only 85 runs 
resulted in an evolutionary branching while 15 runs had no branching during 
the simulation time. Panel c: N = 10 5 , /i = 10~ 8 , out of total 100 total only 81 
runs succeeded and 19 had no branching during the simulation time. 



5.2. Location of the first evolutionary branching 

The trait substitution process considerably slows down as x 
approaches zero, so that it must be considered as a sequence of 
discrete jumps of size e towards zero, rather than a continuous 
process. At levels close to zero a more refined analysis of fitness 
asymptotics is required. 

In the following we derive an approximation for the location 
of the first evolutionary branching. We show that it happens 
close to zero, but typically not at x = 0. More precisely, we 
demonstrate that evolutionary branching becomes possible in 
an interval of width of order e around x = 0. It is thus necessary 
to distinguish between trait values separated by a few mutation 
steps. For simplicity we assume in the following that initially 
x > 0, so that the trait substitution process approaches x — 
from above (as in the simulation results shown in Fig. |3}. The 
problem of determining the location of the first branching event 
is symmetric and corresponding expressions for x < are easily 
found. 

The aim is to find trait values x > where mutations to and 
from x - e are mutually invasive. The corresponding condition 
is: Afj._ £ > 1 and M^ _e) > 1 . Using Eq. ( fl3l l and dropping the 
terms much smaller than e 2 we arrive at two conditions: 



M 
M 



■*1 

x—e 
lx-e) . 



1 « xe + ae /2>0, and 
1 * -(x-e)e + ffe 2 /2>0. 



(22) 



We conclude that the first evolutionary branching is possible for 
< x < (1 + a/2)e, so that populations with trait values x and 
x — e are mutually invasive provided: 



e < x < j*e 



with 



fa/21 . 



(23) 



Here \a/2~\ denotes the smallest integer greater than a/2. (To 
avoid complications arising in the situation when j* = 1 + a/2 
and Mi - 1 = o(e 2 ) for x = j*e, we will assume that a is 
not an even integer.) The fact that the critical trait value j*e 
is larger for larger values of a is very intuitive: as competition 
becomes weaker, initially impossible invasion of (J* - 1 )e into 
7*6 becomes favorable once the benefit of reduced competition 
outweighs the cost of having a "worse" trait value. 

Fig. shows where the first evolutionary branchings oc- 
curred in an ensemble of 300 simulations of the stochastic, 
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Figure 6: Shows coexistence of trait values prior to the first branching, a Evo- 
lution of the integer trait value j(t) as a function of time. Parameters: a = 5, 
e = 10~ 2 , N = 5 X 10 6 , ii = 10~ 9 . The branching region is given by — 3 < j < 3. 
b Dependence of fdf) upon time. Clearly seen is the coexistence of the trait 
values j = 2,3 at f\ = 11/12 and = 1/12. In the case shown, coexistence 
lasts sufficiently long for evolutionary branching to occur. 



individual-based model. The parameter a was chosen to be 
a = 9 which implies that j* = 5. The parameters were cho- 
sen so that the overall substitution rate was the same for all 
runs, 9 = pN = 1CT 3 . All simulations were started at xq = le, 
and were run for the same total time f, given by 6t — 10 5 . In a 
small number of runs evolutionary branching did not occur dur- 
ing this time. Fig. [5] shows where the first branching occurred 
for the remaining runs. First, we see that apart from a small 
number of cases, branching occurs in the region predicted. Sec- 
ond, the larger the population size /V is, the more likely it is that 
evolutionary branching occurs at the boundary of the branching 
region (j* = 5 in this case). To explain this behaviour requires 
a more refined analysis. We return to this question below in 
Section l531 Third, in the limit e — > the condition d23b is con- 
sistent with the standard prediction in this limit (that the first 
branching occurs at x = 0). Fourth, some branchings occurred 
one mutation step earlier than expected. This can be explained 
by the fact that occasionally it can take a long time for a slightly 
advantageous mutant to stabilise, long enough for the next mu- 
tation to initiate the first evolutionary branching. 

5.3. Critical state of coexistence prior to first branching 

In the preceding section we have derived conditions de- 
termining where the first evolutionary branching may occur, 
Eq. ( 1231 ). However, results of our computer simulations of 
the stochastic model described in Section [2] also demonstrate 
that while evolutionary branching may occur (and often does 
occur) when these conditions for evolutionary branching are 
met, branching need not necessarily take place immediately. 



This raises the question how long it takes on average in a fi- 
nite populatio n for evolu t ionary branching to occur (a question 
also raised by iBoettigerl (120101) ). What is the probability that 
evolutionary branching happens when the substitution process 
reaches the boundary of the region where the first evolutionary 
branching becomes possible (Fig. 0? Fig. [6] shows a realisa- 
tion of a computer simulation of this situation for a = 5 (cor- 
responding to j* = 3). In Fig. we see how the trait value 
reaches the critical value j* = 3 from above, and how subpop- 
ulations with trait values j = 3 and j — 2 coexist until evolu- 
tionary branching occurs. The corresponding population sizes 
fj are shown in Fig. |6jj. For evolutionary branching to occur, 
the critical state of coexistence must last sufficiently long for a 
favourable mutation to occur. It is possible, however, that the 
critical state of coexistence may spontaneously disappear by a 
fluctuation. We see in Fig. |6jj that the population-size fluctu- 
ations in the coexistence state are substantially larger than the 
population-size fluctuations in the monomorphic states. 

A preliminary analysis shows that while the fluctuations of 
the monomorphic population sizes are of order N~ 1 ^ 2 , the fluc- 
tuations of the subpopulation sizes in the critical state of coex- 
istence are much larger: of order (Ne 2 ) -1 ' 2 in the limit of large 
values of /V and small values of e. The smaller value e takes, 
the shorter is the average life time of the critical state of co- 
existence, lowering the probability that evolutionary branching 
occurs in this state. Without going into details, we just briefly 
sketch our argument (which remains to be made rigorous). Let 
(fk,gk) be the densities (/V ,ffg' J ) of two populations co- 
existing at time k and having the neighbouring trait values je 
and fe with / = j - 1. The evolution of the vector (fk,gk) can 
be approximated by a stochastic dynamical system 



2ft 



1+fk+agk 
Sk+l ~ T+bf k 



1 



Uk, 
Wk, 



e 



-eV+2-2y) 



b = e -<?(a'+2f) > 



(24) 



where and Wk are independent normal random variables with 

2 yjfk<Jk+agt) an( j 2 y/gt(bft+gt) 



zero means and standard deviations 



i+fk+agk l+bfk+gk 

respectively. Provided the deterministic part has a nontrivial 
stable point we consider a linearised version of this system 
around this stable point. Our analysis indicates that while the 
sum fk + gk behaves as a stationary autoregression process with 
fluctuations of the order N~ l l 2 , the difference /{ - gk behaves 
as a stationary autoregression process with fluctuations of order 
(A^e 2 ) -1 ' 2 . Given that the fluctuations in the sizes of coexist- 
ing populations featuring the trait values (fe, (j* - l)e) are of 
order (Ne 2 Y 1 ^ 2 a rough estimate of the probability of the sud- 
den loss of coexistence is obtained using the approximate tail 
probability for the normal distribution c\e~ ClNe ~ , where c\ and 
C2 are positive constants. Thus, the life time T of the critical 
state of coexistence is expected to be inversely proportional to 
the probability of sudden loss of coexistence c\e~ C2Ne ~ so that 



logT ~ S a Ne 2 



(25) 



for some constant S a depending on a through (123V Fig. [7] 
shows simulation results for the life time T as a function of 
/V and e. The results of the simulations are consistent with the 
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Ne 2 

Figure 7: Shows time 7" cx[ to extinction of the coexistence state at the boundary 
of the branching region for the first branching, in the absence of mutations 
(symbols). Parameters: /i = 0, a = 3, e = 0.01, 0.02, and 0.04. Unless 
otherwise indicated, two-sigma error bars are smaller than the symbol size. 
The error bars were determined from 1000 independent simulations, with three 
exceptions: for e = 0.01 and Ne 2 = 40, 100 independent simulations were 
used, for Ne 2 = 50 and 60, 20 independent simulations. Also shown are fits to 
the law (25), lines. We find S 3 as 0.052. 

expectation Eq. (l25l l. We find 5 3 « 0.052. The precise mathe- 
matical derivation of Eq. d25T l and the calculation of the constant 
S a are interesting questions for further work. 

5.4. Evolution after the first branching 

Condition Eq. ( |23] > demonstrates that if the first evolution- 
ary branching occurs, then it happens within a region of size e 
around x = 0. We now turn to the evolution of the pair of trait 
values (x\,X2) after the first evolutionary branching. Standard 
theory suggests that the dimorphic population evolves symmet- 
rically (along the dashed line in Fig. |4}. Our simulations of 
the stochastic, individual-based population model with discrete 
mutational steps show that the pair of coexisting trait values fol- 
lowed through the chain of consecutive replacements does not 
develop in a fully symmetrical fashion. The random path on 
the ;ti-X2-plane (see Fig.|4]i follows the trajectory of a random 
walk with steps from (x\ , X2) to either (x\ , X2 + e) or (x\ - e, xi) 
depending on which of the two branches the next replacement 
has been successful. The repulsion between the two branches 
is due to the pressure of competition. It is advantageous to stay 
further apart reducing interspecies competition. 

Observe that the consecutive steps in this random walk are 
weakly negatively dependent: the further the walk deviates 
from the diagonal x\ = —x%, the larger is the probability for 
the next step to decrease this deviation. To see this, consider 
a pair of coexisting populations having trait values X\ and xi 
with x\ < < X2 and for definiteness X2 > \x\\. For a given 
branch, the corresponding replacement rate is a product of two 
factors (see the discussion leading to Eq. ( fTBI l in the monomor- 
phic case): the stable population size and the probability of fix- 
ation for a single mutant. We show that given the lower branch 
is closer to zero, X2 > \x\\, both factors of the replacement rate 
for the lower branch are larger making the move of the lower 
branch (always off zero, see Appendix B) more probable. In- 



deed, as the negative trait value is closer to zero, its stable pop- 
ulation size must be larger than the size of the population with 
the positive trait value. On the other hand, as computed in Ap- 
pendix B, the probability of fixation for a single mutant in the 
X\ -population is also larger, Eq. (l49l l. 

The last observation implies that with high probability the 
locations of the pair of branches satisfy 

X! = -x 2 + 0( Vi) (26) 

prior to the second evolutionary branching. This means that the 
the expected location of the second evolutionary branching, if 
any, should not deviate from the diagonal more than a constant 
times Indeed, the predicted deviation from the diagonal can 
only be made larger if we neglect the negative dependence and 
consider a symmetric random walk with independent steps. Re- 
ferring to the classical central limit theorem for the simple sym- 
metric random walk we observe that the quantity X\ +X2 (which 
is always zero in the symmetric case) is of order e -\fn, where n 
is the number of jumps in the random walk (representing the 
number of replacements since the first evolutionary branching). 
Since the second evolutionary branching takes place at values 
of x\ and X2 of order unity, one must wait for a large number of 
jumps, namely of order 1/e for the second branching to occur. 
We conclude that the path taken by the two branches in the x\- 
X2 plane will deviate from the diagonal x\ = -X2 by a distance 
of at most of order e y/l/e = -\[e. 

5.5. 'Triggering cross 'for the second branching 

We now discuss the location of the second evolutionary 
branching. As in the previous section, this location is deter- 
mined by mutual invasibility analysis, but now for the dimor- 
phic case D = {x\, X2\. The second evolutionary branching may 
occur on either of the two branches stemming from the first evo- 
lutionary branching. We first discuss the condition for the case 
where the second evolutionary branching occurs on the upper 
branch of the dimorphic population (an example is shown in 
Fig. [3ji) . With z — X2 ± e the corresponding conditions for mu- 
tual invasibility are M^ 1 '* 2 ' > 1 and A4* > 1. Analysis of the 
signs of Mi* 1 '* 2 ' - 1 and Mjj* 1 ' 2 ' _ 1 ■> reported in Appendix B, 
shows that in view of Eq. ( TZoT i the second evolutionary branch- 
ing of the upper branch becomes possible in the region 

x\ = -x a + 0( Ve). X2-x a + 0( Ve), (27) 

where x a is given by Eq. (fTTI) . More precisely, we show (see 
Appendix C) that the second evolutionary branching of the up- 
per branch with high probability occurs in the region 

(xi + x a )c a - x 2 + x a = 0(e) (28) 

which is a neighbourhood of the straight line 

(xi + x a )c a =x 2 - x a . (29) 

The coefficient c a is given by 

(l+2a) 2 log(l+2a)-2a(l+a) 

c a ■ (30) 

(1 + 2a)(2a - 1) log(l + 2a) + 2a(l + a) 
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Figure 8: Locations of the second evolutionary branching in the x\ -xj -plane. 
Shown are the lines (29), solid blue, and M2\ , solid red, as well as the line cor- 
responding to symmetric evolution of (jti.xi) (black dashed). Results of sim- 
ulations determining the locations of second branchings are shown as points. 
Blue points correspond to cases where the second branching occurred on the 
upper branch, red points to branchings of the lower branch. Parameters: a = 9, 
e= l(r 2 , N = 4x 10 6 ,/i = 2.5x KT 10 ,#i<4x 10 3 . Six sets of initial condi- 
tions were used: 0*1. /z) = (-20, 30), (-30, 20), (-10, 35), (-25, 25), (-35, 10), 
and (-1,45). 



It is shown in Appendix A that c a satisfies c a < 1. The 
minimum value of c a equals 0.7732 . . . and is achieved at 
a = 4.0533.... 

Turning to the possibility that the second branching occurs 
on the lower branch we find a similar condition 



(X2 X(^)C a X\ 



0(6) 



(31) 



corresponding to an e-neighbourhood of another straight line 

X\ + x a = (x 2 - x a )c a . (32) 

Taking these results together, we have shown that the second 
evolutionary branching is expected to occur in the vicinity of 
the lines Eqs. d29l and d32b . If the pair (x\,x?) of diverging 
branches comes close to the line Eq. 1291 . then a second evolu- 
tionary branching of the upper branch becomes possible. Con- 
versely, when the pair (xi , xi) comes close to the line Eq. ( 1321 . a 
second evolutionary branching may occur on the lower branch. 
Note that since // is assumed to be independent of the trait value, 
our answers Eqs. (129) and (l32l depend only on a. The cross 
formed by the pair of lines (l29l . and (l32l is centered at the 
point {-x a , x a ) in the xi-X2-plane. 

Eqs. (|29l , (l32l are consistent with results of direct nu- 
merical simulations of the model. Fig. [8] shows the lines 
d29l and d32l for a — 9 and e = 10~ 2 , as well as the loca- 
tions of second evolutionary branchings for an ensemble of 
simulations of the stochastic, individual-based model. Six 
sets of dimorphic initial conditions were used, (xi,X2) = 
(-0.2, 0.3), (-0.3, 0.2), (-0.25, 0.25), (-0.35, 0.1), (-0.1, 0.35) 
and (0.01, 0.45). We remark that the last three initial conditions 
deviate substantially from the diagonal (dashed line in Fig. @J. 
As shown above, paths in the xi-X2-plane typically exhibit de- 
viations of order -\fe from the diagonal. We have nevertheless 
included dimorphic initial conditions far from the diagonal in 



order to clearly exhibit the locations of the second evolutionary 
branching. 

In Fig. [8] blue dots correspond to cases where the second evo- 
lutionary branching occurred on the upper branch of the dimor- 
phic state, while red dots correspond to cases where the second 
evolutionary branching occurred on the lower branch. The co- 
ordinates of the dots show the location of the second branching 
event in the xi-X2-plane (we determined this location by taking 
the average of the trait values of the critical state of coexistence 
prior to the second branching, as well as the coordinate of the 
second branch when coexistence first occurred). Fig. |8]is con- 
sistent with the theoretical expectation that second branchings 
of the upper branch are triggered by the line Eq. (1291 . while the 
second branchings of the lower branch are triggered by the line 
Eq. (32). 

6. Conclusions 

We have analysed evolutionary branching in a stochastic, 
individual-based model for the evolution of a trait value subject 
to small but non-negligible mutational step sizes e. We found 
that the branching patterns are in general asymmetric at distinct, 
non-null mutational steps , complementing th e standard adap- 
tive dynamics predictions iGeritz et al.l (119981) valid in the limit 
e — » 0. In particular we found conditions describing the loca- 
tions where the first branching may occur, Eq. (l23l . and where 
the second branching may take place, Eqs. ( 1271 . ( l28l . and (OH . 
Results of simulations of a stochastic, individual-based model 
at small mutational step sizes e were seen to be consistent with 
these conditions. Our results are derived in the limit of large 
population sizes N and small mutation rates //, but allow for 
fixed mutational step sizes. In this respect our results Eqs. (l23l , 
( f27t . ( 1281 . and d3TT > complement the established approach. 

We have also demonstrated that population-size fluctuations 
in the critical state of coexistence prior to a branching may erase 
this state. Indeed, its sojourn time depends sensitively on mu- 
tational step size. Evolutionary branching occurs typically only 
when this time is much longer than the time between successful 
mutations. 
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Appendix A 

This appendix is divided into three sections summarising ma- 
terial needed in appendices B and C. 

1. We discuss the properties of the three functions: 



C(a) = --{{\+aY- a l )xi, 

( 1 + 2a) 2 
c i( ff ) = . ~ lo g(! + 2a) - 1, 



2a(l + a) 
(1 +2a)(2a- 1) 



log(l + 2d) + 1 



(33) 
(34) 
(35) 



2a(l + a) 

Here x a is defined by Eq. (fTTT i. We show that C(a) > 0, and < 
c\{a) < c 2 {a) for all positive values of a. The latter inequalities 



in view of c a - c\(a)/c 2 (a) imply that the constant c a , defined 
by Eq. (l30b and determining the slopes of the triggering lines 
Eqs. ([29), d32l l, satisfies < c a < 1. 
To see that C(a) > we observe that 

a l+2a h(l+2a) 



2 4(1+ or) 



8(1 + a) ' 



where h(x) — x 2 - 2x log x—l. The function h(x) takes positive 
values for all x > 1, since h(l) = and h'(x) > for x > 1. 
Thus C(a) > for all positive a. 

Similarly, the fact that < c\(a) < C2{a) follows from 

fti(l+2g) 4 

ci(a) = — — , c 2 (a) = ci(a) + -C(a) , (37) 

2a{ 1 + a) a 

where h\(x) = x 2 log x— is positive for x > 1 (since hi(l) = 
Oand h\(x) > 0). 

2. We discuss properties of the equilibrium densities in the 
dimorphic case D = {x\,X2\- In this case, the equilibrium den- 
sities f£ uX2] and ji* 1,X2] satisfying the system Eq. ©: 

w ( ../r-^+.c- 5 "- 1- 

can be represented as fx* 1 = <f>(x\,X2) and /If 1 '* 2 ' = (f>(x2, x\) 
with 



4>(xi,x 2 ) = 



1 - Rx,XtRx 



I _ e -(xz-xi)(axz-(2+a)xi) 
1 _ e -2(\+a)(x,~x 2 ) 2 



(38) 



In particular, in the symmetric case x\ = -x and x 2 = x, we 
obtain 



4>(—x, x) — <p(x, —x) 



1 



I + e -m+a)x 2 



(39) 



3. We quote the Taylor expansion for the function defined by 
Eq. © (recall that Mf = 1 for x e D) 



M D - 1 



(ax-(l+a)A°)A (40) 
+ (| - a 2 x 2 + 2a(l + a)xA D K - (1 + a) 2 B°^ A 2 
valid for small values of A = z - x. Here the terms 



yeD 



assume the form of first and second moments of the coexist- 
ing trait values using the weights Rxyfy satisfying condition 
Eq. ©. 

Appendix B 

In this appendix we discuss in which region of the x\-x 2 - 
plane the second evolutionary branching becomes possible. It 
may occur on either of the two branches stemming from the 
first evolutionary branching. We start by discussing the evolu- 
tion of the upper branch of the dimorphic population on its way 
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towards the second branching (an example is shown in Fig.|3^). 
With z = x 2 ± e the corresponding conditions for mutual invasi- 
bility are M { * uXl] > 1 and M { Xl uz] > 1. The signs of ifcf^ 1 '* 2 ' - 1 
and M£ - 1 are determined by Eq. (l40b with 



A { £" X2] = X2C/>(x U X 2 ) + XiR X2Xl cf>(X2,X 1 ) 

= xi +(x 2 - xi)<p(x u x 2 ), 
= x 2 + (x 2 - JCi)0(JC1 , x 2 ) , 



(41) 



(42) 



where <p(xi,Xz) is given by Eq. ( TJST l. Using Eqs. (1401 . (14T1 . and 
(l42l we find 



M ta,* 2 ) _ ! ^ Ci(jci, Jta)(z - x 2 ) + C 2 (jci, x 2 )e 2 , (43) 



where 

Ci(x u x 2 ) = ax 2 -(l +a)x x (44) 

- (1 + a){x 2 - xi)cj){xi,x 2 ) , 

C 2 (x u x 2 ) = ^-{ax 2 - (1 + a)xi) 2 (45) 

+ (1 + a){x 2 -xi){{a-\)x 2 

- (1 + a)xi}cf>(xi,x 2 ) . 

Using a counterpart of Eq. ([TBI for the survival probability 
P X2+6 (xi,x 2 ) of a single (x 2 + e)-mutant in a stable (xi,x 2 )- 
dimorphic population we find that 

P x , +e {xi,x 2 ) * 2e{x 2 -xi){a-{x 2 -xiY l xi -(1 + a)(p(x u x 2 )). 

(46) 

Similarly, consider the possibility that the second branching 
occurs on the lower branch. In this case we have for z — x\ + e 

M {x,M _ ! w Ci(j(a,jci)(z- jci) + C 2 (x 2 ,xi)e 2 (47) 

which yields 

P Xl - £ (x u x 2 ) m 2e(x 2 - x\)(a + (x 2 - x l )~ 1 x 2 (48) 
-(1 + a)4>{x 2 ,xi)). 

Close inspection of Eqs. (l46l and ( l48l for xi < < x 2 using 
Eq. ( 1381 reveals that 

P„- e (xi,x 2 ) > P X2+e {xux 2 ) given |xi| < x 2 . (49) 



As explained in Section 15741 this is one of the factors ensuring 
the negative dependence among the steps of the random walk 
in the xi-x 2 -plane leading to the condition Eq. (l26l , 

Applying Eq. (1261 we can further refine our previous analysis 
and deduce from Eq. (l43l 



M 



1 * Ci(jci,*2)(z - x 2 ) + C 2 (-x 2 , W ■ (50) 



Replacing x\ by -x 2 in this expression incurs error of order yfe. 
This implies: 

M {*i.*2) _ J = Cl{ _ X2t X2){z _ X2} + 0(e 3/2) _ (51) 



It follows from Eq. (|39]l that 



Ci(— x, x) = (1 + 2a)x ■ 



2(1 + a)x 

I + e -4(l+ff).t 2 



(52) 



a- ,o 4a( 1 + a)x 2 

C 2 (-x,x) = --(! + 2 ff ) 2 x 2 + 1+ ^ 4(1+ ; g);c2 . (53) 

Now let x„ be given by Eq. (fTTT i. Since Ci(-x, x) is positive for 
x e (0, Xq.), Eq. ( BTT ) confirms that the replacement process on 
the upper branch goes upward until x 2 = x a + 0( -\fe). 

Similarly, consider the possibility that the second branching 
occurs on the lower branch. In this case using Eqs. (l26l and 
d47b we get a lower branch counterpart of Eq. (ISTt 



M u„, 2 ) _ j = Cl (xi,-xi)(xi -z) + Oie 3 ' 2 ) 



(54) 



implying that the replacement on the lower branch goes down- 
ward until X] = —x a + 0( Ve)- We conclude that the second 
evolutionary branching on one of the two branches is only pos- 
sible in the region Eq. (1271 . 

Appendix C 

In this appendix we establish Eq. d28l l in Section 1331 which 
in turn implies Eq. d29l determining the triggering line for the 
upper branch. Given Eq. (l27l . the deviations 6i = x\ + x a and 
6 2 -x 2 - x a satisfy <5,- = 0{ Ve). We demonstrate first that 



1 



1 



Ci(x l5 x 2 ) = -ci{a)5\_- -c 2 {a)5 2 + c n {a)5\ 

+ ci 2 {a)8\5 2 

+ c 22 (a) d\ + o(e) , 



(55) 



where c,(a) are given by (f34-b and (1351 . while the exact form of 
the constants Cij(a) does not matter. Using </>(-x a , x a ) = 2 || 2 ^ 
and the Taylor expansion of Eq. d38l we obtain 

1 + 2a S\ 

0(xi,x 2 ) = — +8\<p\ +6 2 (p 2 + — 0ii 

2(1 + a) 2 

5\8 2 S\ 
+ ^ L 4>\2 + —4>22 + o{e), (56) 

where 0,- = 4-(/>(xi,X2)\ Xl =- Xa ^ 2 = Xa and <pjj are the correspond- 
ing second order derivatives. It follows in view of Eq. 
that Eq. (J56j holds with 

ci(a) = -1 - 4(1 + a)x a <pu 
c 2 (a) = 1 + 4(1 + a)x a (f>2, 
cii(ff) = (1 +ar)(0i -*„0n), 

ci2(ar) = (1 + a)(4>2 ~ <P\ ~ x a <pn), 
c 22 (a) = -(1 + a)(<p 2 + x„^ 22 ), 



(57) 



so that it remains only to verify that the new relations for c,(a) 
agree with Eqs. (l34l and (l35l . For this it is sufficient to note 
that 0i = and 2 = ■ This follows from 

Eq. d38l . We conclude that, given Eq. ( |27| ). we have 

2Ci(xi,x 2 ) = ci(ff)(*i + x a ) 

-c 2 (a)(x 2 - xa) + 0(e) . (58) 
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To conclude our derivation of Eq. (TfFt we notice that in ac- 
cordance with Eqs. (l50l l and (l27l i we have 

m {m,xi) _ l a c 1 (x u x 2 )(z - x 2 ) + C(a)^ , (59) 

where C(a) = C2(—x a , x a ) is a positive constant given by 
Eq. ([33j. Eq. (ED implies that the inequality M^+'f 1 > 1 is 
equivalent to C\(x\,X2) > —eC(a). Similarly, due to Eq. ( |5T1 >. 
the inequality Mi* 1 '* > 1 translates into C\(x\, x 2 ) < eC(a). 
Combining the last two relations with Eq. d58l l we derive 
Eq. ( l28l l. In other words, second evolutionary branching on 
the upper branch of the dimorphic population occurs in an e- 
neighbourhood of the line given by Eq. (29). 

Using Eq. d47i i and repeating the arguments summarised 
above, we find the condition (31) for the second evolutionary 
branching to occur on the lower branch of the dimorphic pop- 
ulation. This branching occurs in an e-neighbourhood of the 
straight line Eq. (32). 
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